function fig3_20

% code to solve the BVP
%  y'' = f(x, y, y')  for  x0 < x < x1   '
%  y(x0) = y0  ,    y(x1) = y1

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.05,'P',0.04)

% set boundary conditions
x0 = 0.0;	y0 = 3;
x1 = 1.0;	y1 = 5;

% parameters for calculation
nx = 4000;
error = 0.000001;

% start off with a linear solution
x=linspace(x0,x1,nx+2);
for ix=1:nx+2
    y(ix)=y0+(y1-y0)*x(ix);
end
dx=x(2)-x(1);
dxx = dx*dx;

nep=400;
epM=0;
epm=-4;
eps=linspace(epm,epM,nep);
for ip=1:nep
    ep=10^eps(nep-ip+1)

    err=1;

    counter=0;
    while err > error

    	% calculate the coefficients of finite difference equation
    	a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
    	for j = 2:nx+1
    		jj=j-1;
    		z = (y(j+1) - y(j-1))/(2*dx);
    		a(jj) = 2 + dxx*fy(x(j), y(j), z,ep);
    		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z,ep);
    		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z,ep);
        end

    	% Newton iteration
    	v(1) = v(1)/a(1);
    	u(1) = - (2 + c(1))/a(1);
    	for j = 2:nx
    		xl = a(j) - c(j)*u(j-1);
    		v(j) = (v(j) - c(j)*v(j-1))/xl;
    		u(j) = - (2 + c(j))/xl;
        end
    	vv = v(nx);
    	y(nx+1) = y(nx+1) + vv;
    	err = abs(vv);
    	for jj = nx:-1:2
    		vv = v(jj-1) - u(jj-1)*vv;
    		err = max(err, abs(vv));
    		y(jj) = y(jj) + vv;
        end
    	counter=counter+1;
    end

    newton_iterations=counter

end

% BL solution
nxx=30;
xx=linspace(0,0.1,nxx);
for ixx=1:nxx
    yy(ixx)=3*exp(-xx(ixx)/sqrt(ep));
end

% plot computed solution
plot(x,y,'-','LineWidth',1.3)
hold on
%plot(xx,yy)
grid on
box on
axis([-0.02 1.02 1 5])

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('Solution','Interpreter','latex','FontSize',14)

hold off

%exportgraphics(gcf,'/Users/mark/Desktop/hw243.pdf')


function g=f(x,y,z,ep)
g=(-tanh(z)+y-1)/ep;

function g=fy(x,y,z,ep)
g=1/ep;

function g=fz(x,y,z,ep)
g=-(sech(z))^2/ep;

































